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Abstract 

This  paper  examines  the  role  of  continuity  of  the  basis  in  the  computation  of  turbulent 
flows.  We  compare  standard  finite  elements  and  NURBS  (non-uniform  rational  B-splines) 
discretizations  that  are  employed  in  Isogeometric  Analysis  [23].  We  make  use  of  quadratic 
discretizations  that  are  C'°-continuous  across  element  boundaries  in  standard  finite  ele¬ 
ments,  and  Cl -continuous  in  the  case  of  NURBS.  The  variational  multiscale  residual-based 
method  [2,  4,  12,  22,  38]  is  employed  as  a  turbulence  modeling  technique.  We  find  that  C1- 
continuous  discretizations  outperform  their  (^-continuous  counteiparts  on  a  per-degree- 
of-freedom  basis.  We  also  find  that  the  effect  of  continuity  is  greater  for  higher  Reynolds 
number  flows. 
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1  Introduction 


Due  to  the  complicated  multiscale  nature  of  turbulence,  numerical  simulation  of 
turbulent  flows  remains  a  significant  computational  challenge  to  this  day.  A  com¬ 
munity  of  researchers  focused  on  understanding  fundamental  physical  behavior  of 
turbulence  through  numerical  simulation  typically  makes  use  of  very  simple  ge¬ 
ometries  and  high-order  spectral  or  compact  finite  difference  methods  (see,  e.g., 
[32,  34]).  The  underlying  function  spaces  utilized  in  spectral  methods  are  of  high 
continuity  ( C°°  in  the  cases  of  Fourier  series  and  global  polynomials).  On  the  other 
hand,  turbulent  flows  are  also  of  great  interest  in  general  geometry  industrial  ap¬ 
plications.  These  are  typically  computed  using  finite  volume  and  finite  element 
methods,  which  employ  low-order  approximation  functions  that  are  at  most  C°- 
continuous. 

Recently,  Hughes  et  al.  [23]  introduced  a  new  computational  technique  termed  Iso¬ 
geometric  Analysis  in  an  attempt  to  generalize  and  improve  on  the  finite  element 
method  in  the  areas  of  geometry  and  solution  representation,  as  well  as  mesh  re¬ 
finement.  Isogeometric  analysis  is  based  on  the  geometric  primitives  of  computer 
graphics  and  computer-aided  design  systems.  The  first  instantiation  of  isogeomet¬ 
ric  analysis  made  use  of  NURBS  (non-uniform  rational  B-splines,  see  [16,  36,  37]). 
NURBS-based  isogeometric  analysis  was  successfully  applied  to  solve  problems 
of  vascular  fluid-structure  interaction  in  [5].  Accuracy  of  the  new  methodology  for 
structural  vibrations  was  shown  to  be  superior  to  that  of  standard  finite  elements  in 
[17].  A  mathematical  theory  of  NURBS-based  isogeometric  analysis  was  put  forth 
in  [3].  A  distinguishing  feature  of  isogeometric  analysis  is  so-called  ^-refinement, 
in  which  the  order  of  the  functions  is  increased  together  with  their  continuity.  As 
a  result,  isogeometric  analysis  allows  for  higher-order  and  higher-continuity  dis¬ 
cretizations  on  complex  geometries, 5  and  may  be  thought  of  as  bridging  the  gap 
existing  between  the  procedures  employed  in  the  flow-physics  and  industrial-flow 
communities. 

A  variational  multiscale,  residual-based  turbulence  modeling  framework  has  re¬ 
cently  emerged  as  a  new  concept  in  large-eddy  simulation  (LES)  of  turbulent  flows 
(see  [2,  4,  12,  22,  38]).  In  this  methodology  one  obtains  an  exact  equation  for 
the  coarse,  or  resolved,  scales  and  their  explicit  dependence  on  the  fine,  or  unre¬ 
solved,  scales,  and  the  modeling  task  amounts  to  accurately  representing  the  fine 
scales  in  the  coarse-scale  equations.  In  this  paper,  inspired  by  the  theory  of  sta¬ 
bilized  methods,  we  employ  simple  algebraic  models  for  the  fine  scales,  which 
involve  appropriately-scaled  residuals  of  the  underlying  partial  differential  equa- 


5  It  is  important  to  note  that  the  highest  possible  order  of  continuity  of  the  solution  space 
in  a  NURBS-based  isogeometric  analysis  is  limited  to  the  continuity  of  the  basis  used 
in  the  definition  of  the  geometrical  domain  of  interest.  Pure  ^-refinement  with  maximal 
smoothness  is  only  attainable  in  simple  geometries 
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tions.  For  a  summary  of  the  early  literature  on  stabilized  methods  see  Brooks  and 
Hughes  [10].  Recent  work  on  stabilized  methods  is  presented  in  [1,  8,  9,  1 1,  15,  18- 
21,28,33,41], 

Spectral  studies  of  simple  advective  and  diffusive  model  problems  in  [4]  indicated 
better  accuracy  for  the  G'1 -continuous  case,  especially  in  the  upper  part  of  the  spec¬ 
trum.  NURBS-based  isogeometric  analysis,  in  conjunction  with  the  multiscale  ap¬ 
proach,  was  applied  to  turbulent  flow  computations  in  [2].  Preliminary  results  in¬ 
dicated  better  accuracy  for  higher-order  and  higher-continuity  discretizations.  It 
was  noted  in  [2]  that  when  going  from  a  C°-continuous  linear  to  a  C1-continuous 
quadratic  B-spline  basis  an  increase  in  solution  accuracy  was  obtained.  The  ques¬ 
tion  this  raises  is  whether  it  was  the  order  or  continuity  of  the  basis  functions  that 
led  to  superior  accuracy.  It  is  the  primary  goal  of  this  paper  to  answer  this  question. 
It  should  be  noted  that  good  accuracy  of  B-spline  discretizations  for  wall-bounded 
turbulent  flows  was  also  noted  in  [29-31,  40]. 

The  variational  multiscale  formulation  employed  in  this  study  utilized  the  “advec¬ 
tive  form”  of  the  convection  term  rather  than  the  integrated  by  parts  “conservative 
form”.  The  advective  form  is  often  employed  in  finite  element  flow  codes,  perhaps 
more  so  than  the  conservative  form.  In  earlier  works  on  residual-based  approaches 
[2,  4,  12,  22],  the  conservative  form  was  employed. 

The  paper  is  organized  as  follows.  In  Section  2,  we  state  the  weak  formulation  of 
the  incompressible  Navier-Stokes  (INS)  equations  at  the  continuous  level.  In  Sec¬ 
tion  3,  we  give  the  semi-discrete,  residual-based  variational  multiscale  formulation 
of  INS  in  advective  form.  In  Section  4,  we  present  numerical  results  for  equilib¬ 
rium  turbulent  channel  flows  at  friction-velocity  Reynolds  numbers  180  and  590. 
We  compare  C°-  and  C'1 -continuous  quadratic  elements  and  assess  accuracy  on  the 
basis  of  the  number  of  degrees-of-freedom.  High-fidelity,  direct  numerical  simula¬ 
tion  (DNS)  results  of  [35]  are  used  as  reference  solutions.  In  Section  5,  we  draw 
conclusions. 


2  Incompressible  Navier-Stokes  Equations 


We  begin  by  considering  a  weak  formulation  of  the  incompressible  Navier-Stokes 
equations.  Let  V  denote  the  trial  solution  and  weighting  function  spaces,  which  are 
assumed  to  be  the  same.  We  also  assume  u  =  0  on  T  and  Jn  p(t )  dQ  =  0  for  all 
t  6  ]0,  T[.  The  variational  formulation  is  stated  as  follows:  Find  U  =  {u,p}  £  V 
such  that  VW  =  {w,  q }  £  V, 


B(W,U )  =  ( W,F ) 


(1) 


3 


where 


B(W,  U) 


{Vw,u®u)n  +  (q,  V-u)Q 
(V-  ™,P)n  +  (Vsw,2vVsu)n, 


(W,F)  =  (w,f)a, 


(2) 

(3) 


and 

Vsu  =  ^(Vu  +  VuT).  (4) 

/  is  the  body  force  (per  unit  volume)  and  v  is  the  kinematic  viscosity. 

Variational  equations  (l)-(3)  imply  weak  satisfaction  of  the  linear  momentum  equa¬ 
tions  and  incompressibility  constraint,  namely 

Qu 

—  +  V  ■  (u  <g>  u)  +  Vp  —  V  ■  (2 uVsu)  —  f  =  0  in  (2,  (5) 

V  ■  u  =  0  in  (6) 


Note  that  one  may  use  the  incompressibility  constraint  to  simplify  the  momentum 
equation  as 

du 

— — b  u  ■  Vw  +  Vp  —  uAu  —  f  =  0  in  fl.  (7) 


3  Multiscale  Residual-based  Formulation  of  INS  Employing  the  Advective 
Form 


We  consider  a  direct-sum  decomposition  of  V  into  “coarse-scale”  and  “fine-scale” 
subspaces,  Vh  and  V',  respectively, 

v  =  vh  ®  V  (8) 

Vh  is  assumed  to  be  a  finite-dimensional  space,  which  will  be  identified  later  with 
the  space  of  functions  with  which  we  actually  compute.  One  obtains  a  unique 
decomposition  in  (8)  with  the  aid  of  a  linear  projection  operator  V,  that  gives 
uh  =  VU  e  Vh  and  U'  =  (I  -  V)U  e  V  from  a  given  U  e  V  (see  [25] 
for  details). 

By  restricting  the  weighting  space  to  Vh  in  (1)  and  employing  the  direct-sum  de¬ 
composition  (8)  for  the  solution  space,  we  obtain  the  equation  system  for  the  large 
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scales,  namely:  Find  Uh  =  {uh,ph}  E  Vh  such  that  MWh  =  { wh ,  qh}  E  Vh, 

B(Wh,Uh  +  U')  =  (Wh,F).  (9) 


(9)  indicates  precisely  the  manner  in  which  the  large  scales  depend  on  U'  = 
WiP'}- 

Combining  (2)  and  (9),  we  obtain:  Find  Uh  E  Vh,  such  that  VWh  E  Vh, 

(™h’d4)n-(Vwh’uh®uh)n+^-U% 

+  (“’’’  IS")  a  ~  ^Wk’ U"  ®  ~ 

n-(v<u'««')o  +  (,‘,v.u')n 

-fV-wh,p')  +  (vswh,2vWsu')  =0  (10) 

V  J  it  V  J  ^ 

For  purposes  of  modeling  the  fine  scales  we  make  the  simplifying  assumption 
(wh,  =  0.  We  note,  however,  that  it  has  been  shown  in  [14]  that  it  is  benefi¬ 

cial  to  incorporate  this  effect  in  modeling  the  fine  scales.  The  term  (\/swh ,  2Vsu'j 
may  be  omitted  by  selecting  a  projector  that  enforces  the  orthogonality  of  the  coarse 
and  fine  scales  in  the  semi-norm  induced  by  this  term  (see,  e.g.,  [2]). 


We  turn  our  attention  to  the  convective  terms  in  (10).  Assuming  incompressibility 
of  the  velocity  field,  namely,  V  •  {uh  +  u')  =  0,  we  compute: 

-  (Vwh,  uh  ®uh)n~  (V«A  uh  <g)  «') 

-  (Vw'\  u'  ®  uh^  -  (\7wh,  u'  ®  u7) 

=  —  (Vwh,  uh  ®  (uh  +  w'))n 

-  (' Vwh,  u'  ( Vmh,  u  <g>  v!) 

=  (wh,(uh  +  u')  ■ 

-  (Vwh,u'  C)w/l)n  -  (Vwh,u'&  u')  (11) 


At  this  point  we  assume  that  is  partitioned  into  a  set  of  subdomains,  such  as  finite 
elements  or  NURBS  elements,  and  on  this  partition  we  have  a  finite  dimensional 
space  of  functions  with  local  support  that  forms  our  approximation  space  for  Uh 
and  Wh.  Let  x  =  denote  the  coordinates  of  element  K  in  physical  space, 

and  let  £  =  (£i}f=1,  denote  the  coordinates  of  element  K  in  parametric  space.  Let 
x  =  £e(£)  :  K  — >  K  be  a  continuously  differentiable  map  with  a  continuously 
differentiable  inverse. 
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We  model  the  fine  scales  as  in  [4]: 


U'  «  -rR(Uh):  (12) 

where  r  is  a  4  x  4  matrix  (in  three  spatial  dimensions)  and  R(Uh)  is  a  4  x  1  vector 
that  collects  momentum  and  continuity  residuals  of  the  Navier-Stokes  equations, 

R(Uh)  =  {rTM(u\ph),rc(uh)}T,  (13) 


in  which 


duh 


rM{u\ph )  =  — -  +  uh  ■  Vuh  +  Vph  -  vAuh  -  f, 
at 

rc(uh)  =  V-uh 


(14) 

(15) 


We  define  r  as  follows: 


T  —  diag(TM,TM,TM,Tc)i 


(16) 


where 


tm 


(-4  +  uh  ■  Guh  +  02G  :  G)-1/2, 

vAr  ’ 


(17) 


tc  =  {g-  rMg )  \ 


with  G  a  second  rank  metric  tensor 

dx  dx' 

and  g  a  vector  obtained  from  the  column  sums  of 


g  =  {di} 

d 

9i  =  E1 

3  = 1 


dx 


>]i- 


(18) 


(19) 


(20) 


The  definition  of  tm  in  (17)  is  inspired  by  the  theory  of  stabilized  methods  for 
advection-diffusion-reaction  systems  (see,  e.g.,  Shakib,  Hughes  and  Johan  [39], 
Hughes  and  Mallet  [24]).  The  definition  of  Tc  comes  from  the  small-scale  Shur 
complement  operator  for  the  pressure  (see  [2]  for  a  details). 

Combining  equations  (10)-(12),  we  obtain  our  discrete  formulation:  Find  Uh  G  Vh, 
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such  that  VWh  e  Vh, 


%-) n  +  K'  (“*  -  r"r“>  '  W%  +  («'■' V  '  “'On 

-  (v  •  wh,  ph)  +  (v*w\  2v\7suh)  -  (wh,  /) 

+  (V'  •  Vwh  +  V</\  rMrM)n  +  (v  ■  wh ,  rGV  ■ 

-  (Vtu\  rMrM  ®  rMrM)n  =  0.  (21) 


4  Numerical  experiments  for  turbulent  channel  flow 

4.  1  Problem  setup 


To  examine  the  effects  of  continuity,  we  compute  turbulent  channel  flows  at  Reynolds 
numbers  i?eT  =  180  and  Rer  =  590,  based  on  the  friction  velocity  and  the  channel 
half  width.  To  assess  the  accuracy  of  the  calculations,  we  compare  with  the  DNS 
of  [35], 

The  problem  setup  is  as  follows.  The  computational  domain  is  a  rectangular  box, 
and  the  flow  is  driven  by  a  constant  pressure  gradient  in  the  stream-wise  direction. 
Periodic  boundary  conditions  are  imposed  in  the  stream-wise  and  span-wise  direc¬ 
tions,  commonly  referred  to  as  homogeneous  directions.  A  no-slip  boundary  con¬ 
dition  is  applied  at  the  walls.  The  no-slip  boundary  condition  is  enforced  strongly, 
that  is,  the  discrete  velocity  is  set  to  zero  at  the  walls.  An  alternative  approach  is  to 
enforce  Dirichlet  boundary  conditions  weakly.  This  is  accomplished  by  appropri¬ 
ately  augmenting  the  semi-discrete  equations  (21)  by  terms  that  enforce  the  no-slip 
condition  weakly  (see  [6,  7]  for  details).  Although  the  latter  approach  was  shown 
to  be  superior  to  the  former,  we  did  not  employ  it  in  the  computations  presented 
herein. 

We  employ  quadratic  finite  elements  that  are  G'°-continuous,  and  quadratic  NURBS 
elements  that  are  C1  -continuous  across  element  interfaces.  The  basis  functions  uti¬ 
lized  in  our  computations  are  shown  in  Figure  1 . 

The  semi-discrete  equations  (21)  are  advanced  in  time  using  the  general izcd-a 
method  [13,  26].  We  use  meshes  that  are  stretched  in  the  wall-normal  direction 
according  to  a  hyperbolic  function  to  cluster  points  near  the  wall.  Moreover,  in  the 
definition  of  tm  (17)  we  set  Cj,  a  positive  constant  that  derives  from  an  element¬ 
wise  inverse  estimate  (see,  e.g.,  [27]),  to  36. 

Numerical  results  for  all  cases  are  reported  in  the  form  of  statistics  of  the  mean 
velocity  and  root-mean-square  velocity  fluctuations.  The  statistics  were  computed 
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(a)  C°-continuous  quadratic  elements 


(b)  G'1 -continuous  quadratic  NURBS 

Fig.  1.  Basis  functions  employed  in  homogeneous  directions.  Basis  functions  are  shown 
in  a  univariate  setting.  Three-dimensional  basis  functions  are  obtained  by  taking  tensor 
products  of  one-dimensional  basis  functions.  The  meshes  are  stretched  in  the  wall-normal 
direction. 

by  sampling  the  velocity  field  at  the  mesh  knots  and  averaging  the  solution  in  time 
as  well  as  in  the  homogeneous  directions.  The  meshes  were  chosen  such  that  the 
number  of  degrees  of  freedom  for  both  quadratic  discretizations  are  approximately 
the  same.  All  computational  results  are  presented  in  non-dimensional  wall  units. 


4.2  Turbulent  channel  flow  at  ReT  =  180 


The  domain  size  is  27 r,  2,  and  4/37T  in  the  stream- wise,  wall-normal,  and  span-wise 
directions,  respectively.  The  corresponding  DNS  computation  was  carried  out  on  a 
domain  of  the  same  size  with  128  x  129  x  128  spectral  functions  in  the  stream- wise, 
wall-normal  and  span- wise  direction,  respectively. 
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(a)  Mean  stream-wise  velocity 


y+ 

(b)  Stream-wise  velocity  fluctuations 

Fig.  2.  Turbulent  channel  flow  at  Rer  =  180  computed  on  a  323  element  mesh.  Comparison 
of  C°-  versus  C1  -continuous  discretizations. 

We  carried  out  computations  employing  C°-  and  C1  -continuous  quadratic  dis¬ 
cretizations  keeping  the  number  of  degrees  of  freedom  nearly  the  same  in  both 
cases.  For  the  C°  case  we  used  a  mesh  of  163  elements,  which  gave  32  x  33  x  32 
basis  functions  in  our  discrete  space,  whereas  for  the  C 1  case  we  employed  a  mesh 
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y+ 

(a)  Wall-normal  velocity  fluctuations 


y+ 

(b)  Span-wise  velocity  fluctuations 


Fig.  3.  Turbulent  channel  flow  at  Rer  =  180  computed  on  a  323  element  mesh.  Comparison 
of  C°-  versus  C1  -continuous  discretizations. 


of  323  elements,  which  led  to  a  discrete  space  comprised  of  32  x  34  x  32  basis  func¬ 
tions.  (The  open  knot  vector  construction  is  responsible  for  the  extra  basis  function 
in  the  wall-normal  direction;  see  Hughes,  Cottrell  and  Bazilevs  [23].) 
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Figure  2  illustrates  that  the  mean  flow  obtained  with  the  C'1 -continuous  discretiza¬ 
tion  is  slightly  more  accurate.  However,  both  discretizations  show  good  agreement 
with  the  DNS  result.  On  the  other  hand,  the  fluctuations  are  significantly  better  in 
the  case  of  Cl  quadratics  (see  Figures  2-3). 


4.3  Turbulent  channel  flow  at  ReT  =  590 


We  also  performed  turbulent  channel  flow  computations  for  Rer  =  590  to  examine 
the  effects  of  the  continuity  as  the  Reynolds  number  of  the  flow  is  increased.  For 
this  simulation  the  domain  size  is  27t,  2,  and  n  in  the  stream-wise,  wall-normal  and 
span-wise  directions,  respectively.  The  corresponding  DNS  used  the  same  domain 
size  with  a  resolution  of  384  x  257  x  384  spectral  functions  in  the  stream-wise, 
wall-normal  and  span-wise  directions. 

As  in  the  ReT  =  180  case,  we  present  results  comparing  quadratic  discretizations 
that  are  C°-  and  C1 -continuous  while  keeping  the  number  of  degrees  of  freedom 
nearly  the  same.  Figures  4-5  show  the  results  obtained  using  a  C°  mesh  of  32  x  33  x 
32  basis  functions  and  a  C 1  mesh  of  32  x  34  x  32  basis  functions.  The  difference 
in  the  number  of  degrees  of  freedom  is,  as  mentioned  previously,  attributable  to  the 
open  knot  vector  construction  for  the  NURBS  basis.  These  meshes  are  considered 
coarse  for  this  simulation,  which  is  manifested  by  the  fact  that  the  mean  stream- 
wise  velocity  is  over-predicted  by  both  discretizations,  the  C 1  discretization  being 
more  accurate  (see  Figure  4).  As  in  the  ReT  =  180  case,  the  velocity  fluctuations 
are  significantly  more  accurate  for  C 1  quadratic  NURBS,  as  shown  in  Figures  4-5. 

We  performed  the  same  comparison  study  on  //-refined  meshes,  resulting  in  64  x 
65  x  64  functions  for  the  C°  basis  and  64  x  66  x  64  functions  for  the  C 1  basis. 
Results  are  presented  in  Figures  6-7.  Both  the  mean  flow  and  fluctuations  are  more 
accurate  in  the  case  of  Cl  NURBS,  for  which  the  mean  steam-wise  velocity  is  very 
close  to  the  DNS  result.  From  the  above  results  one  may  infer  that  solution  accuracy 
depends  more  critically  on  the  continuity  of  the  discretization  for  flows  at  higher 
Reynolds  number. 


5  Conclusions 


In  this  work,  we  presented  a  residual-based,  variational  multiscale  model  of  turbu¬ 
lence  that  is  based  on  the  advective  form  of  the  convection  term,  a  popular  choice 
in  finite  element  flow  codes.  We  compared  turbulent  channel  flow  results  using  C°- 
and  C1 -continuous  quadratic  discretizations.  Using  a  Cl -continuous  quadratic  ba¬ 
sis  yields  more  accurate  mean  flow  and  fluctuating  quantities  than  (7° -continuous 
quadratic  basis  functions.  We  conclude  that  smooth  NURBS  basis  functions  have 
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(a)  Mean  stream-wise  velocity 


(b)  Stream-wise  velocity  fluctuations 

Fig.  4.  Turbulent  channel  flow  at  Rer  =  590  computed  on  a  323  element  mesh.  Comparison 
of  C°-  versus  C1  -continuous  discretizations. 


advantages  over  C°-continuous  finite  elements  in  turbulent  flow  calculations  as  an¬ 
ticipated  in  [4] . 
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(a)  Wall-normal  velocity  fluctuations 


(b)  Span-wise  velocity  fluctuations 

Fig.  5.  Turbulent  channel  flow  at  Rer  =  590  computed  on  a  323  element  mesh.  Comparison 
of  C°-  versus  C1  -continuous  discretizations. 
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(a)  Wall-normal  velocity  fluctuations 


(b)  Span-wise  velocity  fluctuations 

Fig.  7.  Turbulent  channel  flow  at  Rer  =  590  computed  on  a  643  element  mesh.  Comparison 
of  C°-  versus  C1  -continuous  discretizations. 
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